In [1]:
%pylab inline
import matplotlib.pyplot as plt
from ipywidgets import *  # az interaktivitásért felelős csomag

#import collections

Populating the interactive namespace from numpy and matplotlib

Phonon spectrum for diamond systems including 1st and 2nd radial derivatives of the pairpotential between nearest neighboring atoms

Result in Michael P. Marder: Condensed Matter Physics (2nd Edition) is reproduced.

see here on page 351 in Figure 13.7.

In [ ]:
In [2]:
# Az abra kimentesehez az alabbiakat a plt.show()  ele kell tenni!!! 

#savefig('fig_rainbow_p2_1ray.pdf');  # Abra kimentese
#savefig('fig_rainbow_p2_1ray.eps');  # Abra kimentese

# Abra es fontmeretek
xfig_meret= 8  #    12 nagy abrahoz
yfig_meret= 6    #   12 nagy abrahoz
xyticks_meret= 15  #  20 nagy abrahoz
xylabel_meret= 20  #  30 nagy abrahoz
legend_meret= 21   #  30 nagy abrahoz
In [3]:
def matafn(h):
    'Array flatten a matrix list of appropriate dimensions'
    H=hstack(h[0])
    for sor in range(1,len(h)):
        H=vstack([H,hstack(h[sor])])
    return H
In [4]:
def Dm_op(qv, f1, f2):
    
    ## q-points are in units of 2 pi/a, qv

    #unit cell vectors (fcc):
    a1 = array([0,1/2, 1/2])
    a2 = array([1/2, 0, 1/2])
    a3 = array([1/2, 1/2, 0])
    
    
    cell = array([array([0, 0, 0]), array([0, 1/2, 1/2]), 
                  array([1/2,0, 1/2]), array([1/2, 1/2, 0])])
    
    delta = array([array([-1/4, -1/4, -1/4]), array([-1/4, 1/4, 1/4]), 
               array([1/4, -1/4, 1/4]), array([1/4, 1/4, -1/4])])


    Dm11 = array(zeros((3,3)))
    Dm12 = array(zeros((3,3)))
    
    for i in range(4):
        d = delta[i]
        a = cell[i]
        # it seems that a factor 1/2 is missing but the results are ok. 
        Dm11 = Dm11 + (f1-f2)*outer(d,d)/dot(d,d) + f2*eye(3)
        Dm12 = Dm12 + ((f1-f2)*outer(d,d)/dot(d,d)+ f2*eye(3))*exp(1j*2*pi*dot(qv,a)) 
        
    
    Dm22 = Dm11
    Dm21 = array(zeros((3,3)))
    
    for i in range(4):
        d = delta[i]
        a = cell[i]  
        Dm21 = Dm21 + ((f1-f2)*outer(d,d)/dot(d,d)+ f2*eye(3))*exp(-1j*2*pi*dot(qv,a))
    
    Dm = matafn([[Dm11,Dm12],[Dm21,Dm22]])
    
    return  Dm
In [30]:
def rajz_phonon(f1_in,f2_in,Npoints):
    
    ## diamond lattice:

    f1cim = f1_in # f1/M in units of (2pi)**2 THz^2
    f2cim = f2_in   # f2/M in units of (2pi)**2 THz^2
    f1= f1cim*(2*pi)**2
    f2= f2cim*(2*pi)**2

    #Npoints = 100;  ## hany pontbol keszul a gorbe
    #ymax = 3.5;

## k-points are in units of 2 pi/a

## high symmetry points in the Brillouin zone
## FONTOS: a szomszedos spec. k pontokat kell venni 
## (ezek most nem a fenti abranak megfelelo pontok, de szomszedos pontok, 
## es igy jo az abra)
    Gp = array([0,0,0])
    Xp = array([0,1,0])
    Wp = array([1/2,1,0])
    Kp = array([3/4,3/4,0])
    Lp = array([1/2,1/2,1/2])
    Up = array([1/4,1,1/4])


#specpoints = array([Gp,Xp,Up,Kp,Gp,Lp])
    specpoints = array([Gp,Xp,Up,Gp,Lp])

    klist = []
    specpos = [0]
    kk = [0]
    tmp = 0

    for i in (range(len(specpoints)-1)):
              
        sl = specpoints[i+1]-specpoints[i]  
        sl_hossz = sqrt(dot(sl,sl))
        dk = sl/Npoints
        for num in arange(0,Npoints):
           # k vector along the line given by speclines:
            klist.append(specpoints[i] + num*dk)
           # length of the k vector and shifted: 
            tmp +=  sqrt(dot(dk,dk)) 
            kk.append(tmp) 
        specpos.append(tmp)

    klist.append(klist[-1]+dk)   

    enlist = []
    for kv in klist:
        enlist.append(sqrt(abs(eigvalsh(Dm_op(kv, f1, f2))))/2/pi)
    
    enlist= array(enlist)

    eigszam= len(enlist[0,:])

#print(klist)
    
# drawing the figure
    figsize(12,7)

    for ii in range(0,eigszam):
        plot(kk,enlist[:,ii],color='r')

#specNev = ["$U$","$\Gamma$","$L$"]
#specpoints = array([Gp,Xp,Wp,Xp,Up,Kp,Gp,Lp])
    specNev = ["$\Gamma$","$X$","$U$","$\Gamma$","$L$" ]
            
    i = 0
    for xc in specpos:
        axvline(x=xc,color='b')
    #text(xc,-0.5,specNev_map[tuple(specpoints[i])],fontsize=xylabel_meret)
        text(xc,-1.,specNev[i],fontsize=xylabel_meret)
        i=i+1   
    
    ylabel(r'$\omega(q)/(2\pi) \,(\mathrm{THz})$',fontsize=xylabel_meret)

#  kikapcsoljuk a xticks: 
    xticks([], [])
    ticks=arange(0,18,3)
    yticks(ticks)

    xlim(0,specpos[-1])
    ylim(0,15.5)

    cim = "phonon spectrum for Si,  $f_1$ ="+str(f1cim)+"$*(2\pi)^2 \, \mathrm{THz}^2$" + \
     ", $f_2$ ="+str(f2cim)+"$*(2\pi)^2 \, \mathrm{THz}^2$"
    title(cim,fontsize=xylabel_meret)

    grid();
In [31]:
rajz_phonon(78,4,10)

Same result as in Michael P. Marder: Condensed Matter Physics (2nd Edition), Fig. 13.7

In [ ]:
interact(rajz_phonon,
        f1_in=FloatSlider(min=0.,max=200,step=10.0,value=78.0,description='f1='),
        f2_in=FloatSlider(min=0.,max=50,step=1.0,value=4.0,description='f2='),
        Npoints=IntSlider(min=10,max=300,step=10,value=50,description='Npoints='));
In [ ]: